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FOREWORD 


In  the  past  three  or  four  years,  several  methods  (e.g. ,  see  References  1  and  2)  have 
been  devised  and  used  to  remove  ephemerls  errors  present  in  geoid  height  and  vertical 
deflection  data  obtained  from  the  GEOS-3  satellite.  The  methods,  though  varying  to  some 
extent,  all  have  two  common  features.  These  are  that  1)  the  satellite's  ground  track  is 
split  into  ascending  and  descending  tracks  and  that  2)  an  a  priori  functional  form  is  assumed 
for  the  ephemeris  errors  along  each  track. 

There  are  several  shortcomings  in  this  approach.  First,  no  a  priori  estimates  of  the 
ephemeris  errors  along  each  track  are  available .  Lack  of  such  information  tends  to  make 
the  selection  of  the  functional  form  haphazard.  Secondly,  correlations  between  time- 
contiguous  ascending  and  descending  tracks  are  not  accounted  for.  Failure  to  do  so  pro¬ 
duces  an  unrealistic  solution.  Finally,  discrete  jumps  in  the  ephemeris  may  be  produced 
along  solution-area  boundaries.  With  regards  to  GEOS-3,  this  latter  problem  has  never 
been  addressed. 

The  technique  described  herein  overcomes  all  of  these  drawbacks.  It  is  believed  that 
the  use  of  this  technique  will  resolve  the  inconsistencies  observed  in  previously  processed 
GEOS-3  data  and  enhance  present  processing  of  SEASAT  data. 


W.  C.  Palmer 
Captain,  USN 
Commanding  Officer 
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—This  report  presents  a  new  technique  for  correcting  these  errors.  "Hie  main  fea¬ 
tures  of  the  technique  are  that  J^the  satellite's  ground-track  is  treated  as  one  time- 
continuous  track  which  repeatedly  intersects  itself  and  that  £)■  no  a  priori  functional 
form  of  the  ephemeris  errors  is  assumed.  A  conjugate  gradient-projection  algorithm 
is  used  to  find  the  unbiased,  discrete  function  of  minimum  weighted  variation  which 
produces  the  geoid  height  differences  observed  at  the  ground-track  intersections. 
After  inspecting  the  discrete  solution,  a  suitably-chosen  continuous  function  can  be 
fitted  to  the  discrete  data  and  used  to  reduce  the  ephemeris  errors.  Preliminary 
analysis  of  SEASAT  geoid  height  data  has  indicated  that  an  altitude  correction  to 
an  accuracy  of  20  centimeters  rms  is  possible. 

Although  this  new  method  is  presently  being  applied  to  satellite  data,  the  scope 
of  the  technique  is  much  wider  than  this.  The  algorithm  can  be  employed  to  reduce 
any  time-dependent,  low-frequency  error  present  in  network-type  surveys  (e.g. , 
oceanographic  shipboard  and  airborne  surveys).  Future  applications  include  the 
removal  of  the  diurnal  variation  present  in  magnetic  surveys  and  the  removal  of  the 
nonlinear  gravimeter  drift  present  in  gravity  surveys. 
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INTRODUCTION 


The  orbit  of  either  the  GEOS-3  or  SEASAT  satellite  is  such  that,  over  a  long  period 
of  time,  the  satellite's  ground  track  repeatedly  intersects  itself.  At  many  of  these  inter¬ 
sections,  there  is  a  data  point  in  the  form  of  an  observed  geoid  height  difference.  Since 
high-frequency,  temporal,  and  environmental  effects  have  been  previously  removed  from 
the  satellite’s  radar-altimetry  measurements,  these  observed  differences  are  attributed 
to  low-frequency  errors  in  the  satellite's  ephemeris.^  Let  f(t)  be  the  negative  of  the  ver¬ 
tical  component  of  the  ephemeris  error  and  let  60  be  the  geoid  height  difference  observed 
at  the  jth  intersection.  We  have^  ^ 

f(t,  )  -  f(t  1  =  61,  j  =  1,2,.. .  ,N 

J  J  (1) 

where  t,  >  t  v  j  and  t,  >  t,  for  i  >  3. 

j  Pj  i  j 

In  order  to  correct  the  ephemeris  errors,  we  must  find  a  solution  y(t)  which  satisfies 
Equation  (1)  either  exactly  or  to  some  predetermined  noise  level. 

Two  facts  complicate  the  situation.  First,  a  bias  in  f(t)  can  never  be  detected  through 
Equation  (1)  since  it  gets  cancelled  out  in  the  subtraction  process.  Secondly,  even  within 
the  equivalence  class  of  all  continuous  functions  with  the  same  bias,  the  solution  of  Equa¬ 
tion  (1)  is  not  unique.  To  overcome  these  difficulties,  we  seek  an  unbiased,  low-frequency 
solution. 

Keeping  this  in  mind,  there  are  two  ways  of  approaching  the  problem.  We  can  either 
assume  a  parameterized  functional  form  y(t, aq.a^  ...,an)  for  the  solution  and  perform 
some  type  of  optimization  over  the  parameters  or  we  can  first  seek  a  discrete  solution 


the  latter  approach. 


Obviously,  any  long-wavelength  errors  will  produce  these  differences.  For  simplicity, 
throughout  this  report  we  will  classify  all  long-wavelength  errors  as  ephemeris  errors. 

The  subscripts  "1"  and  "p"  are  labels  distinguishing  later  points  in  time  from  previous, 
or  past,  points  in  time. 


STATEMENT  OF  THE  PROBLEM 


We  can  obtain  an  unbiased,  discrete,  low-frequency  solution  of  Equation  (1)  by  solving 

O 

the  following  problem.  Minimize  the  weighted0  variation 


2N-1 

wi(yi+1-yi): 

i  =  i 


■j  'j 


and  the  additional  constraint 


2N 

2  yi  -  0 

i  =  i 


(2) 


with  respect  to  the  discrete  function  yj,  i  =  1,  . . . ,  2N  subject  to  the  constraints 

yl,  "  V  =  6I.  J  =  ,  (3) 


(4) 


where^ 


y^  =  y(tJ,  y-,  =  y(t,  ),  y  =  y(t  )  , 
j  j 


(5) 


and 


*i  ■=  l2  < . <  'zN-l  <  t-'N 


(6) 


Problem  (2)  -  (4)  can  be  simplified  somewhat  by  solving  Equation  (3)  for  y^.  and  sub¬ 
stituting  the  result  into  Equations  (2)  and  (4).  Before  doing  so,  we  define  * 


r  if 

k' t, , 


1  if  3  J  5  P.  =  i 


pj  if  =  1 


i  =  1 ..... 2N 


(?) 


and 


<5  =0  j  =  1.....N 


(8) 


3 

The  weights  Wj  >  0,  i  =  1,  ...,  2N-1  will  be  used  to  account  for  the  uneven  spacing 
between  intersections. 

4 

It  should  be  noted  that  for  every  i  there  exists  a  j  such  that  either  lj  =  i  or  p^  =  i. 

3 
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A  new  equivalent  problem  can  now  be  stated. 


Minimize  the  weighted  variation 

:;N- 1 

J  =  i  2  w.  (y.  -  y,  +  i  . . ,  -  6  .  )2 

2  l  ki+i  1+1  1 


(9) 


with  respect  to  the  discrete  function  yp^,  j  =  1, ....  N  subject  to  the  constraint 


l! 


?  V  +  l2  5,.  =  0 

j=l  J  j=l  J 


(10) 


In  vector  notation  this  becomes  minimize  the  quadratic  function 


f(z)  =  A  ZT/^Z  +  ET2  +  C  .  (]]) 

with  respect  to  z,  subject  to  the  linear  constraint 

h(z)  =  dTz  +  e  =  0  ,  (12) 

where 


A  =  (am,n>  <13) 

m  =  1 . N 

n  =  1, . . .  ,N 


with 


<a  -  w 
m,m  p 


+  w 


Pm1 


+  w,  +  w 


1m1 


m  =  1 , . . .  ,N 


(14) 


Since,  from  this  point  on,  we  will  be  dealing  with  the  vector  z,  it  will  be  essential  to 
remember  that  Zj  -  yp.,  j  =  1 . N. 


4 


(15) 


where 


a  =  - y  w.  for  m  f  n 
tn,n  i 

i  =  : 


w,  -  \  .in 


if  Pn.  =  Pn  '  1 
otherwi se 


(16-1) 


A 

W„ 


if  pm  =  Pn  +  1 
otherwi se 


(16-2) 


*  _  j  wp 

w,  -  i  Hm 

o 


if  pm  =  !n  '  1 
otherwi se 


(16-3) 


W„  = 


w 


if  pm  =  1r,  +  1 


n  otherwise 


(16-4) 


A 

Wr 


Wt 


m 


if  1  =  pn  -  1 

m  n 

otherwise 


(16-5) 


A 


if  1  =  Pn  +  1 

m  n 

otherwise 


(16-6) 


A 

W- 


if  1  =  1_  *  1 

m  n 

otherwi se 


(16-7) 


w8  = 


if  1  =  ln  +  1 

m  n 

otherwi se 


(16-8) 


=  "WD  ,  5n  -  l  ”  WD  6p  + 1 

m  Pm  1  pm  1  Hm 

+  w,  _  jC -5-,  '  6i  -  i)  "  w)  ^61  +  i  51  ' 

1  m  m  mm  m 


m  =  1 , . . .  ,N 


(17) 


5 


c 


(18) 


i 

i 

! 

i 

I 


In  the  above  definitions. 


2N-1 

=  7  2  Wi(6i  +  1  -  6.f 

i  =  i 

=  1  m  =  1,. . .  ,N 
N 

e  =  j  2  8^ 

j=i  j 


6c  =  °»  62N+1  ■  °*  wo  °»  W2N  =  0 


and,  in  Equations  (14) -(17), 


w0  =0  if  1  =  p  +  1 

p  m  m 

r  m 


We  now  turn  our  attention  to  the  solution  of  Problem  (11)  -  (12). 


(19) 


(20) 


(21-1) 


(21-2) 


! 
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COMMENTS  ON  THE  SOLUTION  OF  THE  PROBLEM 


We  are  confronted  with  the  problem  of  minimizing  the  quadratic  functionfj 

f(7)  =  j  zTAz  +  bTz 


(22) 


subject  to  the  linear  constraint 


h(z) 


=  d^z  +  e 


=  0 


(23) 


This  problem  has  a  unique  global  solution  providing  that  the  matrix  A  is  positive  definite 


on  the  subspace  Mj  =  {z:c(Tz  =  o] 
the  Lagrangian 


The  problem  can  be  recast  as  that  of  minimizing 


F  ( z ,  A )  =  f(z)  +  Xh(z) 


(24-1) 


=  \  zTAz  +  bTz  +  X(aT?  +  e)  (24-2) 

with  respect  to  the  vector  z  subject  to  the  constraint  (23),  and  the  solution  can  be  found  by 
solving  the  linear  system 


F?(Z,X)  =  0 


(25-1) 


h(z)  =  0 


(25-2) 


Specifically,  we  must  find  the  vector  z  and  the  scalar  Lagrange  multiplier  X  which  satisfy 


(26) 


This  seems  simple  and  straightforward,  and  is  indeed,  providing  N,  the  number  of 
track  intersections,  is  relatively  small.  Large  segments  of  the  track,  however,  may  con¬ 
tain  thousands  of  intersections,  making  the  direct  solution  of  Equation  (26)  extremely 


6  . 

Here,  the  third  term  c  has  been  dropped  from  f(z)  since  it  has  no  effect  on  the  solution. 

7 

The  positive  definiteness  of  A  is  discussed  in  Appendix  A. 


difficult  if  not  impossible.  In  order  to  have  the  capability  to  process  such  large  segments, 
we  shun  any  direct  method  of  solving  (26)  and  select  an  iterative  method,  the  conjugate 
gradient-projection  algorithm.  In  doing  so,  we  will  be  able  to  take  full  advantage  of  the 
fact  that  the  matrix  A  is  sparse  (see  Appendix  A). 


THE  CONJUGATE  GRADIENT-PROJECTION  ALGORITHM 


Conjugate  gradient  algorithms  are  well  documented  in  the  literature.  References  3 
and  4  provide  particularly  lucid  treatises  of  the  subject.  These  algorithms  differ  from 
ordinary  gradient  methods  in  that  at  each  step,  instead  of  moving  in  the  negative  direction 
of  the  gradient,  movement  is  made  along  a  direction  p^  which  is  Q-orthogonal  to  all  previ¬ 
ous  step  directions  p^  i.e. , 


=  0,  1  =  0,1,...  ,k-l  (27) 

O 

where  Q  is  a  symmetric,  positive  definite  matrix. 

For  unconstrained  problems,  the  direction  p^  is  obtained  from  a  linear  combination 
of  the  present  gradient  g^  and  the  previous  search  direction  p. 

^k  =  "^k  +  Vk-i  (28) 


The  directional  coefficient  is  chosen  to  satisfy  the  conjugacy  requirement  (27). 

For  linearly  constrained  problems,  the  direction  p^  is  obtained  as  above  with  the 
exception  that  the  gradient  is  replaced  by  its  projection  onto  the  subspace  defined  by  the 
constraints— hence,  the  nomenclature  gradient-projection.  This  procedure  enables  satis¬ 
faction  of  the  constraints  at  every  iteration. 

To  introduce  our  particular  algorithm,  let  z  denote  a  nominal  point  satisfying  Equa¬ 
tion  (23),  let  z'  denote  a  varied  point,  and  let  p  denote  a  search  direction  which,  when 
multiplied  by  the  stepsize  a,  leads  from  the  nominal  point  to  the  varied  point.  Then,  by 
definition. 


~z  =  ~z  +  a"p  (29) 

Furthermore,  let  the  present  search  direction  p  be  related  to  the  previous  search  direc¬ 
tion  $  via 

p  ■  -r-ft.*)  *  4  (30) 

Equations  (29)  and  (30)  contain  three  unknown  parameters  -  the  Lagrange  multiplier  X. 
the  directional  coefficient  0.  and  the  stepsize  a.  We  will  determine  these  parameters 
so  that  l)  the  constraint  (23)  is  satisfied  at  every  iteration,  2)  the  present  search  direc¬ 
tion  is  A-orthogonal  to  the  previous  search  direction  and  3)  the  maximum  amount  of  de¬ 
crease  in  the  function  f(z)  is  attained  at  each  iteration. 


8 

Vectors  pj,  i  =0,  1,  ...,  k  satisfying  Equation  (27)  are  said  to  be  conjugate  with  respect 
to  Q. 
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r 


For  the  constraint  (23)  to  be  satisfied,  we  must  have 


h(z’)  =  0 

Since  we  have  assumed  z  satisfies  (23),  this  implies  that 


~  Mz  )  =  £  fh(z‘ )  -  h(z) 

=  0 

(32-1) 

(IT?  +  e)]  =  C 

(32-2) 

4  (?•  -  r)T? . 

0 

(32-3) 

=  ~p T  cT  =  o 

(32-4) 

Substituting  Equation  (30)  into  Equation  (32-4)  produces 

(-F£(z,a)  +  0p)Td  =  0 

(33) 

A  T— ^ 

Providing  that  we  have  forced  p  d  to  be  equal  to  zero  in  the  previous  iteration,  this 
reduces  to 

fI(!,a)3  =  o 

(34-1) 

or 

(Az  +  b  +  a3)T3  =  0 

(34-2) 

We  now  choose  X  so  that  Equation  (34-2)  is  satisfied.  We  have 


A  =  -(A z  +  b)Td/3T3  (35-1) 

=  -fJ(3)3/3T3  (35-2) 

■-*Eya 

k 

K-  i 

To  prevent  the  above  logic  from  breaking  down,  on  the  first  iteration  we  set  (3=0, 


(35-3) 


For  the  present  search  direction  p  to  be  A-orthogonal  to  the  previous  search  direction 
p,  we  must  have 

prAp  =  0  (36) 

Coupled  with  Equation  (30),  this  implies  that 

-F-^(z  ,x)Ap  +  SpTAp  =  0  (07-1) 


or 


P  * 


F^(z,x)Ap 

pTAp 


(37-2 


To  achieve  the  maximum  amount  of  decrease  in  the  function  f(z)  at  each  iteration  we 
first  note  that  with  X  given  by  (35-3), 


f(z')»  F(z'.x) 


(38) 


Furthermore,  since  both  X  and  0  have  been  specified,  F(z‘ ,  \  )  is  a  function  of  a  only. 


F(a)  =  F(z* ,X)  (39-1 

=  F(z  +  ap,X)  (39-2; 

=  Z  +  ap)"^A(z  +  ap) 

+bT(z  +  «p)  +  x  [dT ( z  +  up)+  e]  (39-3; 

=  (^p^Ap)a2  +  [(Az  +  b  +  Xd)Tp ]- 

+  i  zTAz  +  bTz  +  x(dTz  +  e)  (39-4; 

=  ( |  PT Ap )  u2  +  p] ( z  ,  \ )  Pa  +  F  ( 0 )  (39-5; 

Minimizing  Equiation  (39-5)  with  respect  to  a  ,  we  have 

F  (  <)  =  pTApa  +  F-T(2 ,\)p  =  0  (40) 


so  that  the  optimum  stepsize  is  given  by 


(in 


ot  =  -F-I(  2,  •  )p7pTA'p 

Note  that  for  a  minimum  we  must  have 

F  (  <)  =  pTAp  •  0  un 


which  will  always  be  satisfied  if  A  is  positive  definite  on  the  subspaee  M^.  (Hoier  liqua¬ 
tion  (32-4).) 

The  bare  basics  of  the  algorithm  have  now  been  derived.  Mam  additional  properties 
of  the  algorithm  exist.  Perhaps  the  most  important  is  that,  when  applied  to  a  quadratic- 
linear  problem  such  as  (22)  -  (23),  convergence  is  assured  in  at  most  N-q  iterations,  where 
q  is  the  dimension  of  the  constraint  h.  Other  properties,  such  as  additional  orthogona lin 
properties,  will  not  be  discussed  here.  We  will  state,  however,  that  due  to  these  additional 
orthogonality  properties,  Equations  (37-2)  and  (41)  can  be  reduced  to 

s  =  Q ( z ,  ) /Q ( z , • )  d.ii 

and 

JC  =  Q(z,x)/pTAp  (441 

where 

Q(z,A)  =  fI(z,A)F-(z,A)  (45) 

The  symbol  " A  "  refers  to  quantities  computed  in  the  previous  iteration. 
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SUMMARY  OF  THE  ALGORITHM 

Tlu-  conjugate  gradient-projection  algorithm  lias  now  been  fully  specified.  The  matrix 
A  and  the  vector  I)  need  to  be  com )xi toil  only  once.  y  After  this  is  done,  the  algorithm  is: 

Step  1,  Select  a  nominal  point  /.  such  that 


h(  z)  =  d^z  +  e  =  0 

Step  2,  Compute 

f-(z)  =  Az  +  b 

Step  5.  Compute 

M 

■-H  2  fz  d) 

k=  k 

Step  1.  Compute 

Fz'(z,-.)  =  f|(z)  +  Ad 


Step  a.  Compute 

Q(z»  )  =  Fl(z,x)F?(?,x) 


Step  6.  Check  for  convergence.  Is 


(46) 


(47) 


(48) 


(49) 


(50) 


(51) 


where  e  is  some  small,  preselected  number  .  If  (51)  is  satisfied,  stop;  if  not,  proceed 
with  step  7, 

Step  7.  Is  the  iteration  number  equal  to  N-q?  If  so,  stop;  if  not,  proceed  with  step  8. 

Step  8.  Compute  the  directional  coefficient  0.  If  this  is  the  first  iteration,  0  =  0;  other¬ 
wise 
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The  matrix  A  and  the  vector  b  require  6  x  N  storage  locations. 
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6  =  Q(z,\)/Q(z,A) 


(52) 


Step  9.  Determine  the  search  direction 


P  =  -F-(z,A)  +  3p 

(53) 

Step 

>  10.  Compute  the  optimum  stepsize 

a  =  Q(z,\)/pTAp 

(54) 

Step 

)  11,  Compute  the  varied  point 

Z  =  Z  +  ap 

(55) 

Step 

>  12.  Replace  the  old  nominal  values  z,  p,  and  Q(z,  X  )  with  the  new' 

nominal  values 

z'  , 

p,  and  Q(z,  X  ) 

z  =  z' 

(56-1) 

P  =  P 

(56-2) 

q(z,x)  =  Q(z  ,a) 

(56-3) 

and  proceed  to  step  2. 

Once  convergence  has  been  attained,  the  remaining  half  of  the  ephemeris  errors  are 
computed  via  Equation  (3). 


II 


ON  CORRECTING  GEOID  HEIGHT  ERRORS  RELATIVE  TO 
A  REFERENCE  SURFACE 


This  section  deals  with  the  situation  in  which  segments  of  *-he  satellite  track  have  been 
previously  corrected  and  it  is  desired  to  correct  the  remaining  segments  relative  to  these. 

In  such  a  situation,  any  of  the  nodes  yij  or  yp^  falling  in  a  previously  corrected  segment 
are  deemed  to  be  errorless.  Let  Iq  be  the  index  set  of  the  errorless  nodes.  Then 

yl  .  =  0  f0r  ]j  €  [o  (57-1) 

J 

y  =  0  for  p.  :  I  (57-2) 

Since 

yi  -  yD  =6,  j  *  i, —  ,n  (58) 

j  j 

we  must  have 

yp.  =  -61.  for  ]j  e  I0»Pj  i  Ic  <59> 

J  J 

Let  us  now  see  how  these  constraints,  namely  Equations  (57-2)  and  (59),  affect  our 
original  problem.  Refer  to  Equations  (9)  and  (10).  We  want  to  minimize  the  weighted 
variation  J  subject  not  only  to  the  constraint  (10),  but  also  to  the  additional  constraints 
(57-2)  and  (59). 

At  this  point,  the  straightforward  approach  would  be  to  substitute  Equations  (57-2) 
and  (59)  into  Equations  (9)  and  (10).  This  would  create  a  simpler,  equivalent  problem  of 
lower  order.  Unfortunately,  it  would  also  create  a  complicated  indexing  problem. 

Without  this  substitution,  Problem  (11)  -  (12)  becomes  minimize  the  quadratic  function 

f(z)  =  j  zTAz  +  bTz  +  c  (60) 

with  respect  to  z,  subject  to  the  linear  constraint 

h(z)  =  CTz  +  g  =  o’  (61) 

where  z,  b,  c,  and  A  are  defined  as  before  and  where  the  N  x  (mi-l)-matrix  C  and  the  (m+1)- 
vector  g  are  given  by 
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0  0  0 


0  0 


i  y 

2^1 

j=. 


1  .  i  f  1  .  I .  and  p  .  ■>  I 


0  if  p.  I 
Ji 


Here,  m  is  the  number  of  errorless  nodes.  The  first  column  of  the  C  matrix  contains  all 
ones;  the  remaining  columns  each  contain  only  one  nonzero  element  whose  value  is  one. 

Problem  (60)  -  (61)  has  a  unique  global  solution  providing  that  the  matrix  A  is  positive 
definite  on  the  subspace  M2  ~  { z :  C^z  =  0  [  .  The  problem  can  be  recast  as  that  of  mini¬ 
mizing  the  Langrangian10 


F(z,  \)  =  f(z)  +  h(z) 


(64-1) 


=  \  z^Az  +  b 1 z  +  xT(CTi  +  g) 


(64-2) 


with  respect  to  the  vector  z  subject  to  the  constraint  (61),  and  the  solution  can  be  found  by 
solving  the  linear  system 


F-(zYx)  =  0 


(65-1) 


h(z)  =  0 


(65-2) 


Again,  the  third  term  c  has  been  dropped  from  f(z)  since  it  has  no  effect  on  the  solution. 
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for  the  vector  z  and  the  vector  Lagrange  multiplier  X  . 

Referring  back  to  the  two  previous  sections,  we  can  trace  the  modifications  to  the  con¬ 
jugate  gradient-projection  algorithm  which  will  be  required  to  solve  the  system  (65).  We 
have: 

P  =  -Mz>a)  +  Sp  f66) 


instead  of  liquation  (50); 


instead  of  liquation  (32-4); 


CTp  =  0 


instead  of  Equation  (33); 


c  ( - Fv( £ , A )  +  sp)  =  0 


C  ( A z  +  b  +  Cx)  =  0 


instead  of  Equation  (34-2); 


instead  of  Equation  (35-1); 


X  =  -(CTC)_1CT(A?  +  $) 


instead  of  Equation  (49). 


F^(z,x)  =  f-£(z)  +  CX 


nr 

Since  CAC  is  an  (m+l)x(m+l)  -matrix,  the  main  concern  now  is  to  try  to  simplify 
Equation  (70)  by  solving  for  the  inverse  of 


C'C  =  1 


.1  I  o 


1  1  1  ....  1 


10 . 0 

0  ■  ’  • 


•  •  0 
.'0  1 


N  r  u 


u  i  I 
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analytically.  We  need  to  find  7  ,  v,  and  W  such  that 


Expanding  we  have 


yN  +  uTv  =  1 

NvT  +  uTW  =  0T 

Y  U  +  V  =  0 

uvT  +  W  =  I 

Equations  (74-1)  and  (74-3)  imply  that 

=  1 

Y  N  -  m 


'1 

1 

.1 


(73) 


(74-1) 

(74-2) 

(74-3) 

(74-4) 

(75) 

(76) 


while  Equation  (74-4)  implies  that 


W  =  I  -  uv 


T  1 


iN  -  m 


N  -  m  +  1  1 

1  • 

1 


1 

1 

1  N  -  *m  +  1 


(77) 
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Substituting  this  latest  result  into  Equation  (71)  yields 


In  compact  notation 

c  Hfz.(s> +  *.)  'j  *  ‘o  a"d  Pj*‘o 

r_  (2,X)  ~  i  J 

J  (o  if  lj£I0  or  pj£I0 

This  completes  the  modifications  to  the  algorithm. 
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SUMMARY  OF  THE  NEW  ALGORITHM 
After  computing  the  matrix  A  and  the  vector  b,  the  new  steps  are: 
Step  1.  Select  a  nominal  point  z  such  that 

h(z)  =  CTz  +  g  =  0 

Step  2.  Compute 

f-(  z)  =  Az  +  ti 

Step  3.  Compute 


(84) 


(85) 


AI  =  - 


N  -  m 


2 

j 


fz.(z> 

J 


Step  4.  Compute  F^(z,  X)  according  to 


(86) 


-(f  (?)  +  >  )  if  1.  i  !,  and  p  $  I 

J  ‘  ^  (87) 

0  if  1  £  I0  or  p.  c  I 

J  w 

Step  5.  Compute 

Q(Z,1)  =  F-T(z,X)F-(z,A)  (88) 

Step  6.  Check  for  convergence.  Is 

Z,X  )/zTZ  <.  e  (89) 

where  e  is  some  small,  preselected  number  .  If  (89)  is  satisfied,  stop;  if  not,  proceed 
with  step  7. 

Step  7.  Is  the  iteration  number  equal  to  N-q?  If  so,  stop;  if  not,  proceed  with  step  8. 

Step  8.  Compute  the  directional  coefficient  0.  If  this  is  the  first  iteration,  0  =  0;  otherwise 


Fz  (Z,A)  = 

j 


8  =  Q(z,A)/Q(z,x) 


(90) 
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Step  9.  Determine  the  search  direction 


P  =  -F-(z,x)  +  ,:p  (91) 

Step  10.  Compute  the  optimum  stepsize 

«  =  Q( z ,x)/pTAp  (92) 

Step  11,  Compute  the  varied  point 

Z1  =  Z  +  :.p  (93) 

Step  12.  Replace  the  old  nominal  values  z,  p,  and  Q(z,  X  )  with  the  new  nominal  values 
z  ,  p,  and  Q  ( z,  X  ) 


z  =  z  * 

(94-1) 

p  =  p 

(94-2) 

Q ( z , A )  =  0(z,A) 

(94-3) 

and  proceed  to  step  2. 

It  should  be  noted  that  when  IQ  is  the  null  set,  the  new  algorithm  automatically  reduces 
to  the  old  algorithm. 


SELECTION  OF  A  NOMINAL  VECTOR 


In  deriving  our  algorithm,  we  have  assumed  that  the  nominal  point-?  satisfies  Equation 
(84).  This  assumption  necessitates  that  we  start  the  algorithm  with  such  a  vector.  Al¬ 
though  any  such  vector  will  suffice,  we  might  as  well  use  our  degrees  of  freedom  wisely 
and  select  a  nominal  that  not  only  satisfies  (84),  but  also  produces  the  desirable  property 
that  the  initial  discrete  function  j'j  ,  y  ,  j  -  1 N  is  of  minimum  modulus . 

j  Pj 

To  accomplish  this,  we  minimize 


J  = 


N 


j=i 


(95) 


with  respect  to  the  vector  z,  subject  to  Equation  (84).  Equations  (84)  and  (95)  can  be 
rewritten  as 


*  [‘V  V  1 

(96) 

and 

N  N 

2  v *7  2  h.  ‘ 0 

j=l  J  j=i  J 

(97-1) 

V=-61.  for  1jeIo*Pj^Io 

J  J 

(97-2) 

yPj  =  0  Pjel0 

(97-3) 

From  Equations  (97-2)  and  (97-3)  it  can  be  seen  that  we  only  need  to  perform  our  minimiza¬ 
tion  over  those  yp.  for  which  lj  i  I0  and  pj  £  Ig,  subject  to  the  single  constraint  (97-1). 
Taking  the  partial  of  the  Lagrangian 

N  N  N 

L  =  J2  [<V  *  6i/  +  JCp.  J*  l(2  «1.)  |i>81 

j=l  J  J  J  j=l  J  j=l  J 

with  respect  to  y^,  lj  i  IQ,  pj  i  Ig  and  setting  it  equal  to  zero  produces 


r  -  2V  +  '1.  +  '  =  0 
Pj  J  J 


or 


■  * '  ^  +  2"  5’j 


while  Equation  (97-1)  can  be  rearranged  as 


N 


j  =  1 . N 

"rV1' 


j  = 

TV" 


2  (»Pj  *  2"!^  -  0 


(99-1  •) 


(99-21 


(100) 


Substituting  (99-2)  into  (100)  results  in 


I  (yp.  *  I  <v  *  fi .>  = 

j  J  J  j  J  J 


Pj^o 


Vro 


(101) 


Using  Equations  (97-2)  and  (97-3)  we  have 

j  =  di'-joi , 


(102) 


where 


1  V  V 

2  j  \i 


-  £  ' 

2  j  'j 


(1  09) 


Placing  this  result  into  (99-2)  yields 


=  _  ^ 

2  +  N  -  m 


j  =  1 . N 

Vd1 


(1114) 


Thus,  our  optimal  nominal  Zj  '  Yp.,  j  -  1,  N  is  defined  via  Equations  (97-2),  (97-3), 

(103)  and  (104).  J 
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When  Iq  is  empty.  Equation  (104)  reduces  to 


as  one  would  expect. 


NUMERICAL  EXAMPLES 


Test  Problems 

The  conjugate  gradient-projection  algorithm  was  developed  into  a  computer-programmed 
subroutine  written  in  FORTRAN  V^The  algorithm  was  then  applied  to  several  test  problems. 
In  each  case,  the  nominal  function!1  was  selected  via  Equations  (97-2).  (97-3).  (!<»( ()  and 
(104).  Convergence  was  defined  by 


v/QU 


Z,A )/zz  <  10' 


(106) 


and  double  precision  arithmetic  was  used  in  the  calculation  of  the  optimum  stepsize  -  and 
the  directional  coefficient  B.  The  algorithm  was  terminated  whenever  Inequality  (106)  was 
satisfied  or  the  number  of  iterations  reached  N-q.^ 


For  these  problems,  the  ephemeris  errors  were  represented  by  various  sinusoids.  In 
each  case,  the  observations  6j. ,  j  =1,  ..  .  ,N  were  corrupted  by  Gaussian  noise  with  a  sigma 
level  of  O.Ol'/iT  (lr-  of  the  maximum  ephemeris  error)  and  the  average  number  of  observa¬ 
tions  contained  in  the  shortest  wavelength  of  the  sinusoid  was  recorded. 


Three  different  weighting  techniques  in  combination  with  three  different  spacing  schemes 
were  used.  The  weighting  techniques  were: 


(1)  Equal  weighting  in  which  all  of  the  weights  Wj,  i  =  1,  . . . ,  2N-1  are  equal  . 

* 

(2)  Inverse  weighting  in  which  Wj  =t - —  ,  i  =  1,  ....  2N-1 

M+l-ti 

where  £  is  a  scaling  factor  and 

$ 

(3)  Inverse  square  weighting  in  which  wi  =  2  ,  i  =  1 . 2N-1  where  again  £ 

scaling  factor.  '  ^ 


is  a 


The  spacing  schemes  used  are  given  below. 


Theoretically,  the  algorithm  has  to  converge  in  at  most  N-q  iterations.  However,  due 
to  round-off  errors,  especially  in  the  calculation  of  the  A-orthogonal  directions,  it  is 
possible  to  impose  a  convergence  criterion  that  is  tighter  than  the  algorithm  is  capable 
of  achieving. 
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Table  1.  Spacing  Schemes 


Scaled  'lime 

Scheme  1 

Scheme  2 

Scheme  3 

.3 

.1 

.20 

.4 

•  *-* 

.  25 

.5 

.5 

.30 

.8 

•  35 

.7 

.y 

.40 

1.3 

i.  l 

•  15 

1.4 

1.2 

.  50 

1.5 

1.5 

,  55 

1.6 

1.  * 

.  60 

1.7 

i.y 

.  55 

2.3 

2.1 

2.20 

2.4 

2  2 

2.25 

2.5 

2.5 

2.  30 

2.6 

2.8 

2.35 

2.7 

2.0 

2.40 

3.3 

3.  1 

2.45 

3.4 

3.2 

2.50 

3.5 

3. 5 

2.  55 

3.6 

3.  N 

2.60 

3.7 

3.9 

2.  55 

. 

. 

4.20 

• 

• 

4.  25 

• 

• 

4.30 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

The  ephemeris  errors  were  as  follows; 


Problem  1: 

Problem  2: 


Problem  3: 


f  ( t )  =  sin  p-t 


f(  t)  -  i[sin  jt  +  sin  .-t  j 


f(t)  '  sin  r-t 


f(t)  =  sin  (jjyt  -  4) 
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Problem  4: 


Problem  >: 


(Hi) 


f(t)  =  2  l s i n  >>-t  +  sin  t| 


Problem  6: 


f(t)  = 


0  for  t  -  t 
sin  -t  for  t  >  t 


In  each  ease.  25  intersections  were  considered.  The  ordered  >et<  p 


and  i.  ;  j 


remained  unchanged  will) 


|p-:j  =  1,...,25|  =  {3,5,8,1,2,6,4,7,9,12,17 ,!*, 11, 14, 22, 

25,24,23,26,29,31 , 32, 35, 34,3$! 


{ 1  . :  j  =  1 . 25}  =  {10, 13, 15, 16, 19,20,21,27, ?fl, 30, 33,36,37, 

38,40,41,42  ,43,44  ,45,46  ,47,48,49,50!  (ni) 

In  Problems  1-5,  the  set  F  was  empty.  In  Problem  6,  it  was  assumed  that  the  track  had 
previously  been  corrected  from  time  t  =  0  to  t  ~  t£5  and  that  there  was  a  huge  time  gap 
between  t95  and  to  (j  -  These  conditions  were  simulated  by  defining  the  ephemeris  error 
via  liquation  (112),  setting  the  weight  w.25  equal  to  zero,  and  defining  -  ]i:  i  -  1, 


The  results  for  these  problems  are  illustrated  in  Figures  1-1>.  Problem  1  is  associated 
with  Figures  1-3,  Problem  2  with  Figures  4-6,  Problem  3  with  Figures  7-'.),  etc.  In  most 
cases  the  conjugate  gradient-projection  algorithm  required  less  than  the  maximum  number 
of  24  iterations  to  converge.  The  plotted  solutions  represent  the  sum  of  the  unbiased  esti¬ 
mate  obtained  from  the  algorithm  plus  the  bias  in  the  true  ephemeris  error.  Note  that  as 
the  average  number  of  observations  per  shortest  wavelength  increases,  so  does  the  accur¬ 
acy  of  the  solution.  At  50  observations/wavelength,  the  accuracy  is  down  to  the  noise  level 
of  the  system.  In  some  cases,  depending  on  the  spacing  scheme,  this  level  can  be  achieved 
with  25  observations/wavelength.  Even  with  as  few  as  5  observations/wavelength  the 
ephemeris  errors  can  still  be  substantially  reduced.  At  this  point,  however,  relative  to  our 
sampling  rate,  the  ephemeris  errors  begin  to  look  like  high-frequency  errors.  Since  the 
conjugate  gradient-projection  algorithm  produces  a  low-frequency  result,  the  generated 
solution  may  be  of  lower  frequency  than  the  true  one,  as  evidenced  in  Figure  14. 


Not  to  be  overlooked  is  the  fact  that  the  algorithm  is  relatively  insensitive  (0  the  weigh¬ 
ting  scheme  employed.  Furthermore,  this  insensitivity  increases  as  the  number  of  obser¬ 
vations  per  wavelength  increases.  Because  of  its  performance  over  a  range  of  observation 
frequencies,  inverse  weighting  was  deemed  to  be  superior  to  inverse-square  weighting. 


► 


SEASAT  Geoid  Height  Data 


Moving  from  "Alice  in  Wonderland"  into  the  real  world,  an  analysis  of  the  first  31)3 
revolutions  of  SEASAT  was  conducted.  From  this  initial  segment  of  the  .-aiell  ite  s  trek, 
good  geoid  height  information  was  obtained  during  125  of  the  revolutions.  From  this  latter 
group  of  revolutions.  25(1-1  ground-track  intersections  containing  data  were  produced.  The 
conjugate  gradient-project  ion  algorithm  was  then  applied  to  the  resulting  2551  geoid  height 
differences.  Since  time  information  had  lieen  prev  iously  removed  from  the  data,  inverse 
weighting  based  on  revolution  number  and  fractions  thereof  was  emplovcd.  The  algorithm 
required  84 u  iterations  to  converge. 

A  cubic  spline  with  evenlv -spaced  nodes  was  then  fitted  in  the  least-square  sense  to 
the  geoid-height -error  output  of  the  algorithm.  The  spline  was  represented  as  a  linear 
combination  of  spline  basis  functions  and  a  pseudo-inverse  solution  of  coefficients  was 
obtained  via  a  singular  value  decomposition  (SYl)1.  Spurious  oscillations  were  damped  by 
zeroing  out  an  appropriate  number  ot  the  singular  values.  A  transformation  was  performed 
prior  to  entering  the  SYD  so  that,  whenever  the  solution  was  non-unique,  the  "minimum 
curvature  "I-  solution  was  obtained  (See  Reference  51. 

The  results  for  revolutions  145  through  151.  207  through  2l:i.  23  i  through  24o.  and  274 
through  280  are  given  in  Figures  10-22.  respectivelv  .  For  these  segments,  the  nns  error 
of  the  spline  fit  ranges  between  17  and  2o  centimeters. 

In  order  to  evaluate  the  accuracy  of  the  geoid  height  correction,  a  fictitious  geoid  height 
error  was  assumed.  As  a  function  of  the  revolution  number  .  this  fictitious  error  was 
defined  by 


f(”) 


3 


sin  37 r 


s  i  n 


2* 


(115) 


This  particular  function  was  chosen  because  its  frequency  content  is  similar  to  that  observ¬ 
ed  in  Figures  10-22  —  the  highest  frequency  present  being  approximately  1'2  cycles  per 
revolution. 


After  corrupting  Equation  (115)  with  Gaussian  noise  at  a  sigma  level  of  10  centimeters, 
fictitious  geoid  height  differences  were  computed  at  the  ground-track  intersections.  These 
differences  were  then  fed  into  the  conjugate  gradient-projection  algorithm;  as  before,  a 
cubic  spline  was  fitted  to  the  geoid-height-error  output  of  the  algorithm.  A  comparison 
between  this  spline  fit  and  the  fictitious  geoid  height  error  is  given  in  Figures  23-26, 

This  comparison  shows  (for  the  revolution  spans  considered  here)  that  over  those  segments 
of  the  satellite's  track  where  no  geoid  height  data  exists  (over  land),  we  can  obtain  an 
extrapolated1,1  estimate  of  the  altitude  error  to  an  rms-accuracv  of  22  centimeters 
(RMS1).  Over  these  segments  containing  geoid  height  data  (over  the  ocean),  we  can 
obtain  an  interpolated  estimate  of  the  altitude  error  to  an  rms-accuracy  of  15  cent¬ 
imeters  (RMS2). 


nr 

13 


The  I_2  norm  of  the  second  derivative  was  minimized. 
Interpolation  between  data  sets. 
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figure  4.  10  ObservotLons/Wovelength 


axis 


.133 


Plgtre  8.  25  Observation* /Wavelength 
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Figure  13.  5  Observotlons/Wovelength 
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Figure  14.  5  Observations /Wavelength 
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Figure  21.  7  Revolutions  of  SERSRT  Data 
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SUMMARY 


A  new  technique  has  been  developed  for  correcting  satellite  ephemeris  errors  indirectly 
observed  from  radar  altimetry.  This  new  technique  differs  markedly  from  previous  methods 
in  that  a  different  problem  is  posed  and  a  different  type  of  solution  (discrete  as  opposed  to 
parameterized-continuous)  is  obtained.  Specifically,  a  conjugate  gradient-projection  algo¬ 
rithm  is  used  to  find  the  unbiased,  discrete  function  of  minimum  weighted  variation  which 
produces  the  geo  id  height  differences  observed  at  the  satellite's  ground-track  intersections. 

It  is  felt  that  this  new  procedure  represents  a  substantial  improvement  over  previous 
efforts.  First,  it  accurately  models  reality  by  retaining  the  concept  of  one.  time-continu¬ 
ous,  satellite  ground  track  which  repeatedly  intersects  itself.  Secondly,  it  takes  into 
account  correlations  between  time -contiguous  ascending  and  descending  tracks.  Thirdly, 
it  allows  for  an  a  posteriori  selection  of  a  functional  form  for  curve  fitting  purposes. 

Finally,  it  requires  a  minimum  amount  of  computer  storage  due  to  the  sparseness  of  the 
coefficient  matrix  of  the  system  to  be  solved. 

It  has  been  demonstrated  that  the  method  is  accurate.  Based  on  the  analysis  con¬ 
ducted  herein ,  it  appears  that  an  altitude  correction  to  an  accuracy  of  20  centimeters 
rms  is  possible. 

The  real  beauty  of  the  technique,  however,  lies  in  the  fact  that  it  is  a  general  method 
that  can  be  used  to  reduce  any  time-dependent,  low-frequency  error  present  in  network- 
type  surveys  (e.g.,  oceanographic  shipboard  and  airborne  surveys).  In  the  near  future,  the 
method  will  be  employed  to  remove  the  diurnal  variation  present  in  magnetic  surveys  as  well 
as  the  nonlinear  gravimeter  drift  present  in  gravity  surveys. 
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APPENDIX  A 


Properties  of  the  A  Matrix 

The  A  matrix  defined  in  Equations  (13)  -  (16)  is  symmetric  and  sparse.  The  sparse¬ 
ness  is  due  to  the  relationship  between  the  sets  {  pj :  j  =  1,  . . . ,  N  (  and  j  1  j :  j  =  1,  . . . ,  N  }  . 
This  relationship  is  such  that,  in  any  given  row  (or  column)  of  the  A  matrix,  there  exist 
at  most  5  nonzero  elements. 

The  matrix  A  is  also  of  rank  N-l.  This  is  due  to  the  nature  of  and  the  relationships 
between  these  nonzero  elements.  From  Equation  (14)  it  can  be  seen  that  each  diagonal 
element  is  equal  to  a  sum  of  weights  and  is,  therefore,  positive.  Furthermore,  from 
Equations  (15)  -  (16)  it  can  be  seen  that  each  nonzero  off-diagonal  element  is  equal  to  the 
negative  of  a  sum  of  weights  and  is,  therefore,  negative.  Not  so  apparent  is  the  fact  that 
the  sum  of  the  absolute  values  of  the  off-diagonal  elements  in  any  row  (or  column)  is  equal 
to  the  diagonal  element  in  that  row  (or  column),  i.  e. , 

N 

am,m  “  ^  jam,n 
n^m 

Theorem  1.  The  matrix  A  is  positive  semidefinite. 

Proof:  From  Gershgorin's  Theorem,  every  eigenvalue  X  of  A  must  satisfy  at  least 
one  of  the  inequalities 

N 

<  2 

n=i 
n^m 


N 

=  2  a  m  =  1 . N  (A-l) 

n=i  n,m 
nsftn 


Substituting  (A-l)  into  (A-2)  yields 


which  implies  that 


X 


<  a 

~  m,m 


m  =  1 . N 


(A-3) 


0  <  A  <  2a_ 


(A-4) 


Hence,  A  is  positive  semidefinite. 


Theorem  2.  The  matrix  A  is  positive  definite  on  the  subspace  -  |x:  xTd  =  0  } 

~T~ 

Proof:  A  is  positive  semidefinite  (and  thus  can  be  split  into  A  A  where  A  is  the  square 
root  matrix).  Hence,  we  only  need  to  show  that  for  nonzero  x  e  Mj,  Ax  ^  o. 
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Suppose  that  xTd  =  0  and  A  x  =  0.  Then  x  e  N(A),  where  N(A)  is  the  null  space  of 
A.  Now  d  e  N(A)  since  Ad  =  6  and  N(A)  is  of  dimension  1  since  A  is  of  rank  N-l. 
Thus  x  =  ad  for  some  at  0.  This  implies  thatx^d  =  ad^d  =  aN  t  0  which  con¬ 
tradicts  the  fact  that  x  c  Mj_.  Hence,  A  is  positive  definite  on  M^, 

Theorem  3.  The  matrix  A  is  positive  definite  on  the  subspace  M2  =  j  x:  C^x  =  0  | 

Proof:  M2C 
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